function idat = fsb_normalize_baseline_divisive(idat,sandbox)

% FSB - DEV : Trial Baseline Normalization by division
%
% EXAMPLE:
% idat = fsb_normalize_baseline_divisive(idat,sandbox)
%
% INPUT:
% idat:         4-D image data
% sandbox:      sandbox experiment struct
%
% OUTPUT:
% idat:         Baseline corrected idat
%
% CALLED BY:
% FSB.m
%
% NOTES:
% Does sometimes create artificial activations in noise due to datatype 
% issues. Be wary when using!
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License: GNU GPL, no express or implied warranties
% 
%$ Revision 1.0
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,'Dividing trials by own baseline intensity...');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Transform int16 to single for further processing
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

idat = single(idat);
idat(idat==0)=NaN;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Define baseline period and do normalization
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

blength=zeros(1,max(sandbox.intrial(:,3)));

for x = 1:max(sandbox.intrial(:,3));

    waitbar(x/max(sandbox.intrial(:,3)));
    scanind = find (sandbox.intrial(:,3)==x);
    if min(scanind)==0;
        scanind(1)=[];
    end
    idat_red = idat(:,:,:,scanind); % Gives volumes in trials
    bsl_trial = sum(sandbox.hemodynamics(scanind,1:4),2); % Doing alternative baseline calculation, stimulation regressors only
    
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Looks for baseline volumes by looking for values that are below
    % 1/20th of maximum response
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    bsl_max = nanmax(bsl_trial);
    bsl_ind = find(bsl_trial<(bsl_max/20));   
    
    bsl_ind2 = nanmean(idat_red(:,:,:,bsl_ind),4);

    blength(x)= size(bsl_ind,1);

    voxel_corr = 1/bsl_ind2; % correction factor voxel by voxel
    voxel_corr_4D = repmat(voxel_corr,[1 1 1 size(idat_red,4)]);
    idat_red = idat_red.* voxel_corr_4D;
    idat(:,:,:,scanind) = idat_red;
    
    disp(['bsl_ind2 : ' num2str(nanmean(bsl_ind2(:)))]);
    disp(['mean voxel_corr : ' num2str(nanmean(voxel_corr(:)))]);
    disp(['min voxel_corr : ' num2str(min(voxel_corr(:)))]);
    disp(['max voxel_corr : ' num2str(max(voxel_corr(:)))]);
  
end;

disp(['Mean baseline length: ' num2str(nanmean(blength))]);
close(h);

max_idat = nanmax(idat(:));
min_idat = nanmin(idat(:));
mean_idat = nanmean(idat(:));

multfac = single(intmax('int16'))/(max_idat-min_idat);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
% Multiply data to be able to continue processing it with Sandbox
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
idat = idat*multfac+mean_idat;
idat = int16(idat);
end
